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ABSTRACT 

Aims. We aim to present simulated chemical abundance profiles for a variety of important species, with special attention given to 
spin-state chemistry, in order to provide reference results against which present and future models can be compared. 

Methods. We employ gas-phase and gas-grain models to investigate chemical abundances in physical conditions corresponding to 
starless cores. To this end, we have developed new chemical reaction sets for both gas-phase and grain-surface chemistry, including 
the deuterated forms of species with up to six atoms and the spin-state chemistry of light ions and of the species involved in the 
ammonia and water formation networks. The physical model is kept simple in order to facilitate straightforward benchmarking of 
other models against the results of this paper. 

Results. We find that the ortho/para ratios of ammonia and water are similar in both gas-phase and gas-grain models, at late times in 
particular, implying that the ratios are determined by gas-phase processes. Furthermore, the ratios do not exhibit strong dependence 
on core density. We derive late-time ortho/para ratios of ~ 0.5 and ~ 1.6 for ammonia and water, respectively. We find that including 
or excluding deuterium in the calculations has little effect on the abundances of non-deuterated species and on the ortho/para ratios 
of ammonia and water, especially in gas-phase models where deuteration is naturally hindered owing to the presence of abundant 
heavy elements. Although we study a rather narrow temperature range (10 - 20 K), we find strong temperature dependence in, e.g., 
deuteration and nitrogen chemistry. For example, the depletion timescale of ammonia is significantly reduced when the temperature 
is increased from 10 to 20 K; this is because the increase in temperature translates into increased accretion rates, while the very high 
binding energy of ammonia prevents it from being desorbed at 20 K. 
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1. Introduction 

The ortho-to-para ratio of molecular hydrogen, H 2 , plays a large 
role in the development of deuterium chemistry at the high den¬ 
sity and low temperature attributed to starless cores. The deuter¬ 
ation sequence begins with the exothermic reaction between HJ 
and HD: 

HJ+HD<—> H 2 D + +H 2 +232K. (1) 

The H 2 D + ion can then donate its deuteron to other abundant 
species, such as CO and N 2 , yielding DCO + and N 2 D + , re¬ 
spectively. However, spin-state effects complicate the reaction 
scheme. Chemical species with multiple protons, or deuterons, 
can exist in different spin configurations - for example H 2 has 
two distinct spin states where the nuclear spin wavefunction 
is either symmetric (ortho-H 2 ; hereafter oH 2 ) or antisymmet¬ 
ric (para-H 2 ; hereafter pH 2 ). The difference in energy between 
the ground states of these two spin states, ~ 170 K (iHugo et al.l 
12009) . can be a very important energy reservoir at the low tem¬ 
peratures of starless cores. Indeed, reaction <|T]) can proceed rel¬ 
atively easily in the backward direction when the ortho forms of 
both H 2 D + and H 2 are involved. Consequently, the ortho-to-para 
(hereafter o/p) ratio of H 2 _is a n imp ortant parameter controlling 
deuterium chemistry (iFlower et al J2006bl) . In addition to deuter¬ 
ation, spin-state effects play a large role in other areas of starless 
core chemistry as well. For example, the formation chain of am¬ 
monia depends critically on the N + + H 2 reaction, which requires 
the presence of oH 2 to proceed efficiently at low temperatures 
dDislaire et al.ll2012h . 


Owing to its importance, spin-state chemistry is now widely 
adopted in numerical chemical models, and has been ap¬ 
plied to all t he stages in the star for mation process, i.e., dif¬ 
fuse clouds ( AlbgrtssongLaD 2014a ), star l ess/prestellar cores 
dWalmsley etjk[ 2004 Flower et al. 2004 fPagani et all l2009fc 
ISipilaet ah 201 Ol |20 131: hereafter SI3) and protostellar sys¬ 
tems dTaauet et all 1201 3f iTaauet et al.l 120141) . However, spin- 
state chemistry is usually discussed in the context of specific 
problems adopting very different physical models, and compar¬ 
ison between the results of different works can be difficult. In 
this paper, we aim to remedy this by presenting relatively easily 
reproducible reference results pertaining to spin-state chemistry, 
in physical conditions appropriate to starless cores. To this end, 
we have developed a state-of-the-art spin-state chemical model 
that includes the spin states of light hydrogen-containing species 
(H 2 , HJ...), and the spin states of the species taking part in the 
formation and destruction networks of ammonia and water. Our 
model also includes the deuterated forms of species with up to 
six atoms. To extend the usability of our results, we present re¬ 
sults for both gas-phase and gas-grain models. 


The paper is organized as follows. Sect.[2]describes our phys¬ 
ical and chemical models in detail. In Sect. [3] we present the re¬ 
sults of our calculations. In Sect. [4) we discuss our results and in 
Sect. [5] we present our conclusions. Appendices lAl to IE1 include 
complementary discussion on our main results, and also present 
additional modeling results. 
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Table 1. Adopted values of the various physical parameters (see However, there are some exceptions. For photodissociation reac- 
text for definitions of the parameters). tions, the rate coefficient is defined as 


Parameter Value 


Tgas - 7’dust 

10K 

4 

1.3 x 10 17 s- 1 

A v 

10 mag 

fl g 

0.1 pm 

Pg 

3.0 gem -3 

n s 

1.5 x 10 15 cm 7 

Ei/Eb 

0.77 

Ri 

0.01 


Table 2. Binding energies (corresponding to a water ice surface) 
of selected species. 


Species Binding energy [K] 


H 

450 

h 2 

500 

C 

800 

N 

800 

Na 

1000 

0 2 

1000 

CO 

1150 

O 

1390 

NH 

2378 

OH 

2850 

nh 2 

3956 

NH, 

5534 

H 2 0 

5700 


2. Model 

2. 1. Physical parameters 

In order to facilitate straightforward comparison of future mod¬ 
eling works against the results presented in this paper, we con¬ 
sider homogeneous models. That is, we fix the values of all phys¬ 
ical parameters but the density of the medium. The model pa¬ 
rameters (introduced below) along with their assumed values are 
presented in Table|T| In all calculations, we set the gas tempera¬ 
ture T gas equal to the dust temperature 7)| L1S |. Although we restrict 
the temperature to T = 10 K in the main text, we present refer¬ 
ence results also for T — 15 K and T = 20K in AnncndixICl 
(these results are briefly discussed in Sect. 14.5I >. 

The binding energies of the various species, corresponding 
to a w ater ice surface, are mainly taken from Garrod & Herbst 
(120061 see also lSinilall201 2'). Table[2] presents the binding ener¬ 
gies of selected species. We assume that the binding energy of a 
deuterated species is equal to that of the corresponding undeuter- 
ated species. This approach has been prev iously adopted by, e.g., 
ICazaux et al.l (l2010l) . lTaauet et akl (1201 3l) and S13. To cover typ¬ 
ical density values associated with starless cores, we present re¬ 
sults at four different densities ranging from 10 3 to 10 6 cm~ 3 . 

2.2. Chemical code 

We us e the gas-grai n chemical code discussed in ISipila et all 
(l2010l) . lSipili d2012l ) and SI3, which uses the rate equation ap¬ 
proach to calculate chemical evolution in the gas phase and on 
grain surfaces. The reaction rate coefficient is, for the majority of 
gas-phase reactions, defined with the modified Arrhenius equa¬ 
tion as 

k = a (77300 K f exp (-y/T) . (2) 


kphoto = a exp (-y! A v ) , (3) 

where Ay is the visual extinction. For cosmic-ray-induced dis¬ 
sociation, the rate coefficient is defined as 


kcR =yi£, (4) 

where £ is the cosmic ray ionization rate. In Eqs. ([3ji and ((TJi. y\ 
and 72 are efficiency factors. 

The abundance of dust grains is calculated as 


*g = 


D A^mol 
"Z T 9 
Hr 2 na gPg 


(5) 


where p mo \ is the mean molecular weight per H atom (1.4); p s is 
the grain material density; R d is the dust-to-gas mass ratio. For 
gas-phase chemical reactions involving a grain and another re¬ 
actant with different electric charge, the rate coefficient is multi¬ 
plied by the “J-factor” takin g into account increased reaction ef¬ 
ficiency due to polariza tion dDraine & SutinllT987l iPagani et al.l 
120091 ISipila eTaDl2()Toh . 

Gas-phase chemistry is linked with grain-surface chemistry 
through adsorption and desorption processes. The adsorption 
rate coefficient of species i is given by 


- VjScr, 


(6) 


where Uj = -y/8 Tr T’gas/ mn, is the thermal speed of species i (k K 
is the Boltzmann constant and m, is the mass of species i)\ S is 
the sticking coefficient, set to unity for all species; cr = na 2 is 
the grain cross section assuming spherical grains with radius a g . 
The adopted desorption mechanism is cosmic-ray induced des¬ 
orption, with rate coefficient given by 

kf es = fio kJ D O0 K) = fi 0 vo./ exp 70 K] , (7) 

where fio = 3.16 x 10 -19 is an efficiency factor; vqj = 
\l2n s koEb iln 2 mj is the charac t eristic vibration frequency of 
species i dHasegawa et al.lfl992l : lHasegawa & Herbstlll993 1. In 
the above, n s and E\,j stand for the density of binding sites on 
the grain surface and the binding energy of species i on the grain 
surface, respectively. For simplicity, we consider only cosmic- 
ray induced desorption in this work, and neglect alternative des¬ 
orption mechanisms such as photodesorption or reactive desorp¬ 
tion. 

The rate coefficient for a grain-surface reaction between 
species i and j is given by 

k^aK^Rf +Rf ff )/n g (8) 

where a is the branching ratio of the reaction; the efficiency fac¬ 
tor Kjj - exp (-7i a /Tdust) °r unity for exothermic reactions with 
or without activation energy (£ a ), respectively; n g is the number 
density of dust grains (see below). The (thermal) diffusion rate 

^chff j s gj ven by 

Rf = ^exp(—£ d , i /7’ dust ), (9) 

where N s = n s Ana 2 is the number of binding sites on the grain 
and E dti is the diffusion energy of species i. The diffusion energy 
is determined by assuming a constant value for the diffusion-to- 
binding-energy ratio E d /E\,. 
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Table 3. Initial chemical abundances with respect to nu, and the 
adopted initial H 2 o/p ratio. 


Species 

Initial abundance 

h 2 

0.5 

He 

9.00 x 10- 2 

HD 

1.60 x 10 -5 

O 

2.56 x 10~ 4 

C + 

1.20 x 10~ 4 

N 

7.60 x 10~ 5 

S + 

8.00 x 10~ 8 

Si + 

8.00 x 10~ 9 

Na + 

2.00 x 10~ 9 

Mg + 

7.00 x 10~ 9 

Fe + 

3.00 x 10- 9 

P + 

2.00 x 10- 10 

Cl + 

1.00 x 10~ 9 

H 2 (o/p) ini 

1.00 x 10~ 3 


We assume that the gas is initially atomic with the exception 
of hydrogen and deuterium which are locked in H 2 and HD, re¬ 
spectively (Table [3}. We u sed the sa me ass umptions about the 
initial abundances as iSemenov et al. i J2010l) except for HD for 
which we adopted n(HD)/«(H) - 1.6 x 10 ~ 5 based on the D/H 
ratio measured in the local ISM dLinskv et al.1119951) . The initial 
H 2 ortho/para ratio is (arbitrarily) set to 1 X 10 5 . 

Finally we note that we do not consider quantum tunneling 
on grain surfaces in the main body of the text. However, refer¬ 
ence results including tunneling are presented in Appendix [E] 
where we also discuss the modifications to the gas-grain model 
required when tunneling is included. 

2.3. Gas-phase and grain-surface reaction sets 

The chemical reaction sets for gas-phase and grain-surface reac¬ 
tions used in this work are reworked versions of the networks 
presented in SI3; the present gas-phase network is based on 
osu_01_2QQ£O] instead of the modified version of osu_®3_20®8 
(iSemenov et alJl20T(ih adopted in SI3. In the present model, the 
spin-state chemistry description of S13 is expanded by adding 
the spin states of the species involved in the formation of water 
and ammonia. Also, the present model contains deuterated forms 
of species with up to 6 atoms, including observationally impor¬ 
tant species such as water, ammonia and methanol. We discuss 
the various updates and additions in detail in the following. 

2.3.1. The spin-state separation routine 

In SI3, a routine was developed to automatically spin-state sep¬ 
arate any reaction involving light hydrogen-bearing species (H 2 , 
H/ and H |). In practice, the routine creates new reactions based 
on pre-determined branching ratios. In S13, the branchin g ra- 
tios for most reactions were deduced using the method of lOkal 
(12004 . in which representations of the nuclear spin wavefunc- 
tions of the various species are used to derive selection rules for 
reactive collisions. The resulting branching ratios correspond to 
pure nuclear spin statistical weights under the assumption that 
the nuclei are completely mixed in the reaction, and thus they 
are likely to be most applicable to highly exothermic reactions 
(l0kall2004)l . For example, the reaction 

H 3 + O OH + + H 2 , (10) 

1 See http://www.physics.ohio-state.edu/~eric/ 


with rate coefficient k, separates into three reactions when the 
spin states of the species containing multiple protons are consid¬ 
ered explicitly: 

H^(o) + O OH + + H 2 (o) 

H+(p) + O—>OH + + H 2 (o) (11) 

H 3 (p) + o-4oh + + h 2 (p). 

In SI3, the branching ratio matrices for some reacting systems 
were derived by hand and applied in the o/p separation routine. 
However, for more complicated reactions involving many pro¬ 
tons and/or multiple products, custom separation rules were ap¬ 
plied (see Appendix A in SI3). The motivation for doing so was 
threefold. Firstly, Oka’s method may not be applicable when the 
exothermicity of a given reaction is very low. Secondly, we did 
not want to follow the spin states of species heavier than H 3 
(with the exception of its deuterated forms), so in any reaction 
where such species (e.g., H 2 S) were present, we assumed that the 
species is in its para form. Thirdly, the branching ratio matrices 
were deduced by hand, which quickly becomes tedious when 
multiple protons and product species are present, even though 
the separation method itself is straightforward. 

For this paper, we have developed a new version of the sepa¬ 
ration routine. In the new version, the branching ratios are calcu¬ 
lated automatically using Oka’s method so that we easily obtain 
the branching ratio matrix of any reaction. Consequently, every 
species with multiple protons is separated into its ortho and para 
(and meta when applicable) states. However, in order to maintain 
relative simplicity of the reaction network, we only keep the spin 
states of the light hydrogen-bearing species, and of those species 
included in the formation networks of water and ammonia (see 
Sect. 12.3.21) . 

When spin species other than those mentioned above are 
present in a reaction, we recombine the reactions over the “un¬ 
wanted” spin species. For example, consider the reaction 

H 3 + CH 2 —> CH 3 + H 2 . (12) 

When separated using spin selection rules, reaction O 
branches into the following set of reactions: 


0 H 3 + oCH 2 


0 H 3 + pCH 2 


pH 4 + oCH 2 


Hi 1 
60 





60 


0 CH 3 + oH 2 
oCHj + pH 2 
pCH 4 + oH 2 
pCH 4 + pH 2 
0 CH 3 + oH 2 
0 CH 3 + pH 2 
pCH| + oH 2 
0 CH 3 + oH 2 
0 CH 3 + pH 2 
pCH 4 + oH 2 
pCH 4 +pH 2 
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iti 

pH+ + pCH 2 —^ 0 CH 3 + oH 2 

5*1 

—^ pCHj + oH 2 

f*i 

—» pCH++pH?. 

To determine the branching ratios of the oH J + CH 2 —» CH) + 
(o/p)H 2 system, i.e., neglecting the spin states of CH 2 and CH|, 
we sum over the oH 2 and pH? production pathways and average 
the rate coefficients on the number of spin modifications of CH? 
(two). A similar treatment of the pH^ + CH? —> CH ) + (o/p)H? 
system finally results in four reactions: 

i*i 

0 H 3 + CH 2 CH++ 0 H 2 (13) 

ik\ 

—>CH^+pH 2 (14) 

13 ^ 

pH+ + CH? CH+ + oH 2 (15) 

23*1 

—* CH 3 + pH 2 . (16) 

The branching ratios in Eqs. (IT3l) to (IT6l > are based on the as¬ 
sumption that the ortho and para forms of CH? are equally abun¬ 
dant. Using the statistical ortho/para CH? ratio (3:1), the four 
branching ratios are 33/40, 7/40, 27/40 and 13/40, respectively, 
while assuming that CH? is completely in the para state yields 
the branching ratios 3/4, 1/4, 3/5 and 2/5. If we assume that the 
reaction proceeds through a proton hop process instead of com¬ 
plete scrambling, we obtain the branching ratios 1 , 0 , 1/2 and 1/2 
(lQkall2004l Table 8 b). Evidently, there are differences in the out¬ 
come depending on the assumed ortho/para ratio of CH? or the 
assumed reaction mechanism (which is not known for the vast 
majority of reactions). We note, however, that heavier molecules 
have a negligible effect on the H? ortho/para ratio for example, 
and our assumptions seem justified in the present case where we 
only consider explicitly the spin states of a limited set of species. 

ft should be noted that the Oka method allows the forma¬ 
tion of oH? from pH?, simply because this is statistically pos¬ 
sible. However, the formation of oH? should be unlikely at low 
temperatures where the collisional energies are small compared 
with the energy difference A E/k = 170 K between the ground 
states of oH? and pH?, unless the reaction is sufficiently exother¬ 
mic. In S13, we assumed that the required exothermicity is not 
reached in general, and hence the network favored pH? forma¬ 
tion. In this paper, we slightly relax the assumptions of SI3, and 
apply the separation routine to all reactions other than charge- 
transfer reactions, where we assume that spin states are con¬ 
served. However, we add an activation energy of 170 K to the 
y coefficient in those reactions where oH? is created by reac¬ 
tants whose spin states are not explicit ly considered (for exam - 
ple, CH 3 +CH? —> C 2 H+ +H?; see also lAlbertsson et al.l2014at) . 
The new approach modifies the results of S13 only very slightly; 
see Sect. 13.11 With these changes, the resultant network should 
be bet ter appl icable at higher temper atures (lAlbertsson et alJ 
120 14bt see also lAlbertsson et al.l2014at) . The combination of nu¬ 
clear spin statistical weights and activation energies for probably 
endoergic reactions approximates the method Hugo et all (120091) 
used to calculate the thermal rate coefficients for the H 3 + H? 
isotopic system, but we note tha_t_a detailed stat e-to-state analy¬ 
sis such as that performed by iHugo et akl d2009t) is required to 
study endothermic reactions consistently. 

We apply the same spin-state separation routine to grain- 
surface reactions. No rate coefficient corrections have been made 
to the surface reactions because in the base surface reaction set 


adopted here (iSemenov et alil20l0 l(. the non-addition reactions 
that create H* (the asterisk denotes a surface species), for exam¬ 
ple H* + H?CO ! —> HCO* + H), already present high activation 
barriers. Of course, oH) can be created by (cr-induced) photodis¬ 
sociations; we assume that these events are sufficiently energetic 
so that oH) and pH) are produced in the statistical ratio 3:1. 

We note that the conservation of nuclear spin may not hold 
for grain-surfac e reactions as there is experimental evidence 
to the contrary dFushitani & Momosell2002 ; lHama & Watanabel 
120131) . We discuss this issue briefly in Sect. 14.41 


2.3.2. Spin-state chemistry of ammonia and water 

In this work, we have included the spin-state chemistry 
of nitrogen-containing species involved in the ammonia for- 
mation chain, rec ently stu died using the Oka method by 
Le Gal et al.l (120 14t see also Flower et al.l l2006al iDislaire et al.l 


20121 : iFaureet al.ll2Q 1 3b fRist et al.ll2013l) . We have checked that 


our o/p separation routine produces th e same branching ratios as 
given in Tables B.l. and B.2. in iLe Gal et al.l d2014l) (see also 
Sect.l4~2l) 

Because of the high exothermicities of the hydrogen addi¬ 
tion reactions on grain surfaces creating NH) and NH 3 (~ 4.2 eV 
and 4.7 eV, re spectiv ely, as calculated from Eqs. (3) and (4) in 
I Allen & Rqbinso nlll977l and using data presented in Appendix C 
of lDu et alj|2012l) . we assume that the high-temperature statisti¬ 
cal branching ratios for the formation reactions. 


NH* + H* 


oNH) + H* 


pNH) + H* 


3/4 


1/4 


2/3 


1/3 


oNH) 


pNH) 

onh; 

pnh; 

pnh;, 


are valid. These branching ratios are, of course, given directly by 
the Oka method. 

We have also included the spin-state chemistry of the species 
involved in the water formation network. The branching ratios 
for the most important reactions in the water network are pre¬ 
sented in AppendixlBl Similarly to the nitrogen chemistry, we 
adopt statistical branching for water formation on grain surfaces: 


3/4 

OH* + H* —> oH?0* 

1/4 

pH?0 


(17) 


2.3.3. Deuteration 

In SI3, a deuteration routine was applied to the OSU network, 
and species with up to four atoms were deuterated. In the present 
work, the deuteration routine is extended to handle species with 
up to six atoms, so that the deuterated forms of important species 
such as methanol and ammonia (whose formation depends on 
NHJ) can be in cluded. The deute r ation routine is based on 
that presented by iRodgers & Millarl ( 1996 ). in which deuterons 
are substituted in place of protons in the various reactions and 
branching ratios are calculated assuming complete scrambling 
(see also SI3). 

The calculation of the nuclear spin-state branching ratios for 
multiply-deuterated species is more complicated than for un- 
deuterated species. For species with three or more deuterium 
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nuclei Oka’s method, based on angular momentum algebra, is 
insufficient because there is no one-to-one correspondence be¬ 
tween angular momentum and symmetry representations. The 
statistical branching ratios between the nuclear spin states of 
species with up to five hydrogen or de uterium nuclei can be 
obtained from Tables 111 and IV of IHugo et al.l (I2009I) . In the 
present paper, the spin-state chemistry in reactions between 
species containing H and/or D only is included using these 
rules, and they are mostly available in the literature. Specifically, 
we include the spin- state che mistry of the HJ + H 2 react¬ 
ing system from IHugo et ahl (120091) . complemented by data 
for other reactions concerning light hydrogen-and-deuter ium- 
bearing sp ecies (H/\ FF>D + , HD etc.) from IWalmslev et alJ 
(120041) and iFlower et alJ (120041) . The rate coefficients for disso¬ 
ciative recombination between H!_and its deuterated forms and 
electrons are taken from lPagani et al.l d2009t) . 

The nuclear spin chemistry of multiply-deuterated molecules 
containing heavy elements is not considered in the present 
model, with the exception that it is necessary to predict the ra¬ 
tio of 0 D 2 and pD 2 released in reactions involving heavy species 
to achieve a properly closed reaction network. Furthermore, 0 D 2 
and pD 2 produced by these reactions can be relevant to the chem¬ 
istry of light molecules containing only H and/or D which are 
treated using the appropriate spin selection rules (see below). 

The deuteration procedure used here directly propagates the 
hydrogen spin-state chemistry into the deuterium chemistry, for 
molecules with up to two deuterium nuclei. For hydrogenated 
species like CH 2 , we assumed that the ortho and para forms 
of CH 2 are equally abundant (see Sect. 2.3.1). However, copy¬ 
ing the reaction scheme from hydrogenated to doubly deuterated 
molecules is equivalent to following the correct nuclear spin se¬ 
lection rules for spin 1 systems with the assumption that 0 D 2 and 
pD 2 are instead present in the ratio 2:1 (i.e., the statistical value) 
in complexes with heavy nuclei, e.g., CD 2 . 

As an example of the spin-state separation of doubly- 
deuterated species, consider the reaction 


D + + CH 2 


D + + CHD 


D + + CD 2 


|feexp(-170/r) 


5*2 

§*2 

5*2 

5*2 exp(-86/T) 



|*2 exp(-86/r) 


CD + + oH 2 
CD + + pH 2 
CH + + HD 
CH + + oD 2 
CH + + pD 2 
CD + + HD 
CD + + oD 2 
CD + + pD 2 . 


Reactions yielding pH 2 or oD 2 are assumed to have no activa¬ 
tion energies, whereas those producing pD 2 are assumed to have 
an activation energy of 86 K, corresponding to the energ y sepa¬ 
rati on between the ground states of pD 2 and oD 2 ( iHugo et al.l 
12009 ). As noted above, the branching ratios for the reactions 
producing D 2 are correct under the implicit assumption that 
[oCD 2 ] : [pCD 2 ] =2:1. However, as in the case of H 2 , reac¬ 
tions involving heavy species have a negligible effect on the D 2 
o/p ratio. Note that for the H + + CHD and D + + CH 2 reactions 
above, one would expect H 2 formation in the o/p ratio 3:1, but 
in the example the ratio is instead 2:1 because the deuterated 
reactions adopt the branching ratios from the undeuterated par¬ 
ent reactions, in this case reactions ( IT9t , where the three-proton 
system leads to an H 2 o/p ratio of 2:1. 

We have verified through testing that the spin separation of 
doubly-deuterated (heavy) species is of little consequence to the 
results presented later in this paper, and therefore the separation 
method used here seems reasonable for the present work. We 
will discuss the spin separation of multiply-deuterated molecules 
in depth in a future dedicated paper. 


H + + CH 2 —U CH + + H 2 . (18) 

Since our program performs the spin-state separation of hydro¬ 
gen species prior to deuteration, the above reaction is first di¬ 
vided into the following two branches: 

H + + CH 2 —> CH + + oH 2 , k = exp(-170/r), 

—> CH + +pH 2 , k= l -k 2 . (19) 

The branching ratio (2/3 oH 2 , 1/3 pH 2 ) follows from the nuclear 
spin statistics under the implicit assumption that pCH 2 and oCH 2 
have equal abundances, i.e. [pCH 2 ] = [oCH 2 ]. Deuterating both 
branches leads in our routine to the following reactions: 

5*2 exp(— 170 /r) 

- CD + oH 2 

CD + + pH 2 
CH + + HD 
CH + + oD 2 
CH + + pD 2 
CD + + HD 


2.3.4. Other updates with respect to SI 3 

In addition to the updates to the o/p separation and deutera¬ 
tion routines described above, the rate coefficients of several 
reactions have been updated. First of all, we now adopt the 
osu_01_2®09 reaction set instead of (modified) osu_®3_20®8, 
so any updates of the base reaction set, including anion chem¬ 
istry, are included in our new_model. The modifications made 
to osu_®3_2®88 bv ISemenov et al.l (12010 ). already included in 
SI3, have been incorporated to the new model as well. 

We have also included new rate coefficient data from the lit- 
eratur e. We adopt the r ate co efficients given in Tables 12 and 
13 in lAlbertsson et al.l (i2013l) - except for reactions involving 
hydrogen and/or deuterium only for which we use other data 
as detailed above. We have also included selected data from 
iLe Gal et alJ (120141), which we summar ize i n TablelB .2| We re- 
fer the reader to Albertsson et al.l (1201 3l) and lLe Gal et al.l (1201 4l) 
for the original data references. 

The impact of the rate coefficient updates is analyzed in 
Sect.im 

3. Results 

In this section, we present benchmarking results for spin-state 
chemistry for different values of density, using gas-phase and 
gas-grain models. All of the results presented below correspond 
to T = 10 K. Results for T = 15KandT = 20 K are presented in 


H + CHD 


H + + CD 2 


|*2 

5*2 

5*2 

1*2 exp(-86/r) 
1*2 
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Fig. 1. Upper panels: abundances of selected non-deuterated (left), deuterated (middle) and grain-surface (both non-deuterated and 
deuterated, right) species as functions of time, as calculated with model T1 (solid lines) or model T2 (dashed lines). Middle panels: 
as the upper row, but comparing models T1 (solid lines) and T3 (dashed lines). Lower panels: as the other rows, but comparing 
models T1 (solid lines) and FI (dashed lines). In all calculations, the density of the medium is set to ii[\ = 10 5 cm -3 . The abundances 
of HrO, NH 3 and NFbD represent sums over their respective ortho and para states. 


Appendix[Cj while results without deuterium and including tun¬ 
neling are presented in Appendices [D] and [E] respectively. The 
data for the various figures presented below are available from 
the authors upon request. 

3.1. Comparison of the new model against S13 

In our new model, there are two main sources of updates: 1) 
our new spin-state separation and deuteration routines and 2 ) 
the rate coefficient updates for various species from the litera¬ 
ture and from the transition from osu_®3_2®08 to osu_®l_2®®9. 
Qualitatively, the new spin-state separation routine may affect 
somewhat the abundances of the various species because of the 
increased emphasis on 0 H 2 in the new separation routine (albeit 
with a 170 K barrier in most cases). However, the extension of 
deuteration from four to six atoms is not expected to have a large 


impact on the overall chemistry. In the following, we study the 
effect of the various updates, by successively introducing new 
updates into the S13 chemical model. The models introduced 
below are also described in Table [4] 

First, we constructed a test reaction set by spin-state separat¬ 
ing the S13 base reaction set (modified osu_®3_28®8) with the 
new routine presented here, but using the old S13 deuteration 
routine and including no new rate coefficient data (model T2). 
The upper row in Fig.Q]plots the abundances of selected species 
as functions of time calculated using model T2 (dashed lines) 
and using the old model of S13 (solid lines; model Tl). In these 
calculations, the density of the medium is set to «h = 10 5 cm~ 3 , 
and the other physical parameters correspond to the values given 
in TableQ] In the right panel, the abundances of the respective or¬ 
tho and para states of H 2 O*, NFQ and NH 2 D* have been summed 
over. It is observed that the new spin-state separation routine 
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t [yr] t [yr] t [yr] t [yr] 


Fig. 2. Abundances of selected species calculated with model F2 (i.e., surface reactions and adsorption/desorption excluded) as 
functions of time at different densities, labeled in the upper left panel. The labels indicate densities of the medium; for example, 10 3 
stands for «h = 10 3 cm“ 3 . 


modifies only slightly the abundances of non-deuterated species 
at early times both in the gas phase and on grain surfaces, while 
the abundances of deuterated species are virtually identical in 
both cases. Hence the new oHi-creating reactions with 170 K 
barriers included here (as opposed to SI3) are of little conse¬ 
quence to the chemistry as a whole (although at high tempera¬ 
tures the situation might be different). 

The middle row in Fig. Q] plots the abundances of selected 
species calculated with a test reaction set constructed by ap¬ 
plying the new spin-state separation and deuteration routines 
to the modified osu_03_2008 reaction set; as before, no new 
rate coefficient data is included (model T3). It is observed that, 
overall, the extension of the deuteration scheme has little effect 
on the chemistry. However, deuteration is enhanced in the new 


model at long timescales because of less efficient HD depletion. 
Inspection of the reaction rates reveals that this effect is caused 
by the dissociation of HDO*, the main grain-surface deuterium 
carrier in our model, whose dissociation produced only OD*+H* 
in S13 owing to a bug in the deuteration routine, but here pro¬ 
duces also OH* + D*. Consequently, more D* is available to form 
HD*, replenishing the gas-phase HD abundance. 

In the bottom row of Fig. [I] we plot again the abundances of 
selected species, now calculated with the full model presented in 
this paper, i.e., including the rate coefficient updates discussed 
in Sect. l2.3.4l (model FI). Evidently, the rate coefficient updates 
lead to significant differences with respect to the previous cases, 
both in the gas phase and on grain surfaces. For example, the CO 
formation timescale is much shorter with the updated rate coef- 
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Fig.3. As Fig. [2] but calculated with model FI. Note that the abundance scales are in some cases different than in Fig. [2] 


7 


7 


7 


7 


Table 4. The models discussed in this work. 


Model 

Description 

T1 

the SiDMe^aL_£2013) model 

T2 

osu_®3_20®8 + new spin-state separation routine; 
old deuteration routine; no rate coefficient updates 

T3 

osu_03_20®8 + new spin-state separation routine + 
new deuteration routine; no rate coefficient updates 

FI 

osu_03_20®8 + new spin-state separation routine + 
new deuteration routine + rate coefficient updates 
(full model) 

F2 

model FI without gas-grain interaction (i.e., no 
adsorption/desorption or surface reactions); 

H-> formation from Kona et ah (20131 


ficients, which is mainly due to numerous changes in electron 
DR rate coefficients. Accordingly, the electron abundance drops 


initially very rapidly in the new model. From an observational 
point of view, the most striking difference between the models 
is the abundance of N 2 FR, which is over an order of magnitude 
lower at late times in the new model, resulting mainly from in¬ 
creased (factor 2-3) electron DR rate coefficients and also from 
changes in neutral-neutral chemistry, with respect to the OSU 
values. The changes is N 2 FD and HCO + abundances propagate 
to their deuterated counterpart ions N 2 D + and DCO + . Also, the 
abundance of ammonia decreases by a factor of 2-3 at late times 
in the new model with respect to SI3. 

Surface chemistry is also affected by the rate coefficient up¬ 
dates. In the new model, atomic hydrogen is at early times pro¬ 
duced mainly by neutral-neutral reactions between light hydro¬ 
carbons instead of electron DR reactions of various ions (owing 
to decreased electron abundance), and there is more atomic hy¬ 
drogen available in the gas phase. The hydrogen is then adsorbed 
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onto grain surfaces, leading to efficient production of, e.g., water 
and ammonia which form mainly by hydrogenation. We stress 
that these effects are caused by changes in the gas-phase chem¬ 
istry, as none of the surface reaction rates have been updated. 
These results emphasize the importance of constraining the re¬ 
action rate coefficients, particularly for key reactions. 

Our chemical model predicts gas-phase and grain-surface 
abundances in physical conditions attributed to starless cores, 
and our model results could in principle be used to interpret 
ice observations. However, direct observations of grain-mantle 
species toward starless cores are extremely hard owing to the 
very large extinctions. Thus, we would need to rely on, e.g., ob¬ 
servations toward embedded Class 0 sources, assuming that the 
ice feature originates in the source envelope. A detailed investi¬ 
gation of modeled ice abundances versus observations is left for 
future work, and we focus our attention on gas-phase species in 
the present paper. 

In what follows, all results have been produced with the 
new ortho/para separation and deuteration routines, and include 
the rate coefficient updates from the literature as described in 
Sect. ITTfl 



3.2. Gas-phase chemistry 

Figure [2] present the results of gas-phase modeling at different 
densities. Here, we have not considered adsorption and/or des¬ 
orption of any species, and the formation of fF, HD and Dt on 
grain surfaces is included in the form of gas-phase reactions as 
in the model of lKong et al.l (120131) . 

Ion abundances decrease significantly with increasing den¬ 
sity, which is a consequence of the increased electron DR rates at 
high density, and of the high abundances of heavy neutral species 
owing to the absence of depletion. The abundances of heavy ions 
follow largely that of H) which is primarily destroyed by re¬ 
actions between abundant neutral species like CO and N 2 , and 
electrons. One example of this tendency is HCO + , which is the 
most important reaction partner of free electrons at late times 
both at low and high densities. This ion is formed in a reaction 
between Ht and CO and destroyed mainly in electron recom¬ 
bination which returns CO. At high densities, the enhanced pro¬ 
duction of free electrons by the cosmic ray ionization of H 2 leads 
to an increase in the electron DR rates (at constant temperature), 
and the HCO + abundance decreases accordingly. N 2 H + is locked 
in a similar cycle (H+ + N 2 —> N 2 H + ; N 2 H + + e~ —> N 2 + H), 
dependent on the abundance of HJ. Also the N 2 H + abundance 
decreases clearly when the density increases. 

Deuteration is also suppressed in the absence of depletion, 
because H I reacts preferentially with the abundant heavy species 
rather than with HD. 

3.3. Gas-grain chemistry 

Figure [3] presents the results of gas-grain modeling at differ¬ 
ent densities, adopting the full model (FI) described above. 
Evidently, the inclusion of depletion significantly decreases the 
abundances of carbon and oxygen-containing species with re¬ 
spect to the gas-phase model. The H^ abundance is much higher 
at late times in the full model with respect to the gas-phase model 
because of the depletion of its main reaction partners (e.g., 
CO) onto grain surfaces. The chemistry of nitrogen-containing 
species is different than that of, e.g., carbon and oxygen because 
nitrogen chemistry depends on slow neutral-neutral reactions. 
On the other hand, the binding energy of (atomic and molecu¬ 


Fig. 4. o/p ratios of NH 3 (upperpanels) and H 2 O (lower panels) 
as functions of time at different densities, labeled in the upper 
left panel. The left panels correspond to the gas-phase model 
(F2), while the right panels correspond to the gas-grain model 
(FI). 

lar) nitrogen is relatively low (Table [2j, so nitrogen chemistry 
can still occur efficiently after other heavy elements have been 
depleted onto grain surfaces. Consequently, an appreciable abun¬ 
dance of N? , the precursor molecule of both NH 3 and N 2 H + (see, 
e.g jFontani et al.ll2012l) . can be present in the gas at late times. 
This is demonstrated by two effects. Firstly, the abundance of 
NH 3 is less than an order of magnitude lower in the full model 
than in the gas-phase model at high density, even though deple¬ 
tion tends to decrease its abundance. Secondly, the abundance 
of N 2 H + , dependent on H j, is orders of magnitude higher with 
respect to the gas-phase model even at high density (we assume 
that ions do not deplete onto grain surfaces, but can recombine 
with negatively charged grains). 

Deuteration increases strongly when gas-grain interaction is 
taken into account. However, a local peak in deuteration is ob¬ 
served around the time that HD starts deplete, in line with the re¬ 
sults of S13. HD depletion is somewhat less severe in the present 
model than in S13 because we do not consider quantum tunnel¬ 
ing. In agreement with previous results in the literature, our re¬ 
sults clearly show that depletion is needed to produce observable 
amounts of, e.g., oH 2 D + and N 2 H + at high density. 

Figures |2] and [2 indicate that the o/p ratios of ammonia and 
water depend only slightly on density, both in the gas-phase and 
gas-grain models. We discuss this issue further below. 

4. Discussion 

4.1. The o/p ratios of water and ammonia 

Figure [4] plots the o/p ratios of NH 3 and H 2 O at different densi¬ 
ties according to gas-phase and gas-grain models. The o/p ratio 
of ammonia is very similar in both models up to n n = 10 s cm" 3 . 
The ratio is larger at ~ 10 4 - 10 5 years at n H = 10 6 cm -3 in the 
gas-phase model, in which one of the main production pathways 
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forNH+atf ~ 10 5 yristheH 3 O + +NH 3 —> NH+ + H 2 0 system 
which forms all three spin variants of NH j. However, ntNH J for 
exam ple only dissociates to oNH 3 (see Table 2 in [Le Gal et al.l 
l2014l l. The detailed behavior of the spin-state chemistry is dif¬ 
ferent in the gas-grain model, where the reaction pathways with 
H 3 0 + are practically absent at t > 10 5 yr because of water de¬ 
pletion. Nevertheless our results indicate that, at late times, o/p- 
NH 3 is ~ 0.4-0.5 regardless of the considered model or assumed 
density. 

The o/p ratio of H 2 0 behaves similarly to that of NH 3 in 
the sense that the largest differences between the gas-phase and 
gas-grain models are observed at «h = K/’cm 3 . The o/p ratio 
of HiO is tied to that of H 3 0 + , which evolves differently in the 
gas-grain model owing to multiple spin-state details (H 3 0 + and 
H 2 C) are connected mainly through H 2 C) 1 + H 2 —» H 3 0 + + H 
and the H, + H 2 0 —> H?0' + H 2 system). The o/p-HiO ratios at 
different densities and in different models are, again, very close 
to each other at late times where we derive o/p-H 2 0 ~ 1.5 - 1.7. 

Unlike the o/p ratio of Hi which is strongly influenced by 
grain-surface chemistry, the o/p ratios of H 2 0 and NH 3 are de¬ 
termined nearly completely by gas-phase chemistry because of 
their high binding energies (see Tablc[Tji. The fact that the o/p ra¬ 
tios of the respective species are very similar in both gas-phase 
and gas-grain models supports this interpretation (Fig. [4]). The 
late-time values of the o/p ratios can be justified by a rate coeffi¬ 
cient analysis of the main formation and destruction paths of am¬ 
monia and water; such an analysis is presented in Appendix lAl 

We note that the late-time H 2 0 o/p ratios derived here 
are much lower than th ose typically observed i n the ISM 
( van Dishoeck et alll2013b see also the discussion in lKeto et alJ 
120141) . The disagreement is probably due to the fact that obser¬ 
vations toward dark clouds with intermediate to high density are 
missing. To investigate this issue, we ran test calculations with 
our gas-grain model at conditions simulating translucent clouds 
(«h = 10 2 cm ~ 3 and T - 10K) with two values of Ay, lOmag 
and 1 mag. We found that at Ay = 10 mag, the H 2 0 o/p ratio is 
~ 1.6, but at Ay = 1 mag the ratio is ~ 3, suggesting that pho¬ 
tochemistry plays a large role in determining the o/p ratio in re¬ 
gions with low visual extinction. This result underlines the need 
for accurate physical models when interpreting observations. 

Am monia is a useful p r obe of the gas tem perature in dark 
clouds (iTafalla et al.l 12004 lJuvela et~aki 120121) . and simulated 
NH 3 line emission profiles can be used to constrain the ki¬ 
netic temperature. However, a proper comparison of simulated 
line emission against a detection of, e.g., pNH 3 requires an 
estimate of the NH 3 o/p ratio. Given that the NH 3 o/p ratio 
changes very little with density in our models, we deduce that 
o/p-NH 3 ~ 0.4 - 0.5 is a good conservative estimate of the ra- 
tio in dark clouds. T his result is consis tent with observation s by 
iPersson et al.l (12012l) and the models of Le Gal et al.l (1201 4 ) (see 
also Sect. E3. 

Finally, we stress that we have not carried out a complete 
parameter-space study of the ammonia and water o/p ratios, and 
that such a study should be taken on in the future. 

4.2. Comparison aaainst fTe Gal et a/J (2014) 

Recently. iLe Gal et al.l ( 2014 1) presented gas-phase modeling re¬ 
sults of nitrogen chemistry in dark clouds. They included in 
their model the spin-state chemistry pertaining to the ammonia 
formation network, as we do here. We have run s ome test c al- 
culations to compare our results against those of ILe Gal et al.l 
(120141) . adopting the same physical conditions and initial chem¬ 
ical abundances for He, C, N and Fe (we choose [C]/[0] = 0.8 


Table 5. Abundances and abund ance ratios predicted by our gas- 
phase model (F2) and that of Le Gal et al.l ( 2014 1. 


Parameter 

Our model 

Le Gal et al. (2014) 

[NH] 

5.9 x 10- y 

1.1 x 10~ 8 

[NHi] 

9.5 x 10-'° 

3.5 x 10- 9 

[NH 3 ] 

3.2 x 10~ 8 

1.4 x 10~ 8 

[NHi] / [NH] 

0.16 

0.3 

[NH 3 ] / [NH] 

5.4 

1.3 

[ 0 NH 1 ] / [pNH 2 ] 

1.8 

2.5 

[oNH 3 ] / [pNH 3 ] 

0.4 

0.6 


and [S] = 8 x 10~ 8 1. [Le Gal et all (120141) do not explicitly state 
the initial H 2 o/p ratio; we assumed a value of 3. 

Table [5] summarizes the results of these tests. In our model, 
there are many cosmic-ray-induced and photodissociation reac¬ 
tions with long timescales, which prevent the system from reach¬ 
ing a true steady-state; the values given in the table correspond 
to late-time chemical evolution (f = 2.0 X 10 7 yr), after which 
temporal changes in chemical abundances are generally small. 
The Le Gal et al. steady-state values have been read off their 
Figs. 4 and 5. Evidently, our model gives lower NH and NH 2 
abundances and a higher NH 3 abundance than the Le Gal et al. 
model, but the ortho/para ratios of both NH 2 and NH 3 are similar 
in both models. 

We also tested an initial [C]/[0] ratio of 0.3, and found that 
in this case the late-time NH and NH 2 abundances are lower by 
one and tyro orders of magn itude, respectively, than those pre¬ 
dicted bv lLe Gal et al.l (120141) . However, the NH 3 a bundance and 
the NHi and NH 3 o/p ratios are again similar to the lLe Gal et al . 
d 20 li) model. 

Our test results indicate that there are minor differences 
in fractional abundances between our model and that of 
iLe Gal et alJ (120141) . The two models predict similar o/p ratios 
for NH 2 and NH 3 , which is expected since we adopt th e sam e 
spin-state chemical description for these species as lLe Gal et al.i 

J20I4I) . 

4.3. Hi formation in gas-phase models 

In the gas-phase results presented above, the formation of Hi 
(and its isotopologs HD and P 2 ) is i ncluded in the form of gas- 
phase reactions (see lKong et al.ll^CHjI) . From Fig.[2]it is evident 
that this approach leads to, e.g., an increase of the Hi o/p ratio 
at late times. Another approach to including the formation of 
Hi and its isotopologs in a gas-phase model is to invoke grain- 
surface formation processes, but letting only atomic H and D to 
be adsorbed onto the grain surfaces. 

We have run test calculations at T — 10 K comparing the two 
approaches. The results of these calculations are presented in 
Fig-E where the abundances of selected species as functions of 
time are plotted at different dens ities, us ing different approaches 
to Hi formation. Evidently, the iKong et akl (12013l) approach is 
practically equivalent to the grain-surface formation approach 
when “normal” binding energies (450 K for H and D) are used 
(the dotted and dash-dotted lines are practically superimposed). 
In the gas-grain model with very low H binding energy (100K), 
H desorbs very fast, leading to inefficient Hi production on the 
grain surfaces regardless of the density. Increasing the binding 
energy increases the average time spent by H atoms on the grains 
and thus leads to more efficient Hi production. 

The choice of the Hi formation efficiency has little effect 
on the abundance of H 2 D + because deuterium chemistry is sup- 
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Fig. 5. Abundances of selected species as functions of time at 
different densities (labeled in the figure), in a gas-phase model 
at T = 10 K, using different approaches to H 2 formation. Dash- 
dotted lines correspond to the gas-phase reactions of lKong et al.l 
(120131) . Solid lines correspond to grain-surface H 2 formation as¬ 
suming only H and D adsorption and with the binding ener¬ 
gies of H and H 2 and their deuterated isotopologs set to 100 K. 
Dotted lines (superimposed on the dash-dotted lines) correspond 
to grain-surface H 2 formation assuming only H and D adsorp¬ 
tion and with the binding energies adopted in the main body of 
this work. The abundances of ammonia, water and H 2 D repre¬ 
sent sums over the abundances of their respective ortho and para 
states. 


pressed owing to the lack of depletion, and on the abundance of 
water because its formation pathway proceeds through exother¬ 
mic reactions. However, the abundance of ammonia is slightly 
affected at late times because its formation starts with the N + +H 2 
reaction which is strongly endothermic when pH 2 is involved 
dLe Gal et al-EoIl . 


We have checked that the o/p ratio of ammonia is hardly 
affected by the choice of the H 2 formation efficiency, despite 
its effect on the total NH 3 abundance. This is because the o/p 
ratio is mainly determined by reactions further up in the am¬ 
monia formation chain (see Appendix^. The o/p ratio of wa¬ 
ter is not affected because it is not sensitive to the H 2 o/p 
ratio. However, the o/p ratio of H 2 D + is affected because it 


is mainly determined by the H? o/p ratio ( Paganietak 1 9921 


iGerlich et al .1120021 iHueo et afll2007t iBriinkenet al.l 1201 4l) : the 

late-time H 2 D + o/p ratio changes by a factor of < 2 depending 
on the choice of binding energies. 


4.4. Conservation of spin in grain-surface reactions 

In this work, we have assumed that high-temperature statisti¬ 
cal branching ratios are applicable to both gas-phase and grain- 
surface reactions. However, it has been previously suggested that 
different selection rules ma y apply for grain-surface reactions 
dFushitani & Momosell2002l) and, in the case of H 2 formation. 



Fig. 6. Abundances of selected species as functions of time in a 
gas-grain model at T - 10 K and n h = 10 s cm -3 . The solid lines 
correspond to model FI, while the error bars correspond to test 
models in which grain-surface reactions can only create ortho 
species (upper limits) or para species (lower limits). 


that H 2 ortho/para co nversion may occur on the surface subse¬ 
quent to H 2 formation (IWatanabe et al.l2010UHama & Watanabel 

i2on . 

To constrain the possible effect of uncertainties in grain- 
surface spin chemistry on our results, we constructed two test 
models in which grain-surface reactions can create either only 
para, or only ortho states - with the exception of H 2 and D 2 for 
which we assume the high-temperature statistical formation ra¬ 
tios (i.e., both ortho and para are created). The result of this test 
is shown in Fig. [ 6 ] Evidently, the effect of reactions other than 
the fundamental H* + H* reaction (and of its deuterated coun¬ 
terparts) on the gas-phase spin-state chemistry is negligible, as 
both the upper and lower limits correspond almost exactly to our 
normal full model. 

We have carried out this test for all the medium densities 
considered in this paper («h = 10 3 to 10 6 cm -3 ; with T — 10 K), 
and we conclude that the possible grain-surface branching ra¬ 
tio uncertainties related to spin-state chemistry do no t modify 
the results of this paper. However, we note that I Watanabe et al.l 
(120101) discussed H 2 ortho/para conversion subsequent to H 2 for¬ 
mation; such conversion processes might translate to modifica¬ 
tions in the gas-phase chemistry as the H 2 molecules desorb. We 
currently have no means to test this, i.e., no data to form a rate 
coefficient. This issue should be investigated in the future. 

4.5. Results at higher temperatures 

In Appendix O we present modeling results analogous to those 
shown in Figs. [2] to [4] but calculated assuming either 7„ as = 
T dus t = 15KorT gas = Tdust = 20 K. No dramatic differences be¬ 
tween the temperatures are seen in the gas-phase models, reflect¬ 
ing the weak temperature dependence of ion-molecule chem¬ 
istry. 

However, large differences between the temperatures are 
present in the gas-grain models. For example, nitrogen chem¬ 
istry is much faster at high density at T — 20 K than at T = 10 K. 
The higher temperature increases the adsorption rates, but on 
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the other hand the very high binding energy of ammonia for 
example means that practically no desorption occurs even at 
T — 20 K. Consequently, ammonia and N 2 I-C deplete almost 
totally from the gas phase in a timescale of ~ 10 6 years at 
«h = 10 6 cm -3 . Depletion occurs also at 7’ = 10 K, but in a much 
longer timescale (Fig. [3}- We note that in the case of ISFFD, the 
“depletion” is caused by the freeze-out of N 2 onto grain surfaces, 
as we assume that molecular ions (e.g., N 2 FD) dissociatively re¬ 
combine with negatively charged grains and the products (in this 
case N 2 + H) are immediately returned to the gas phase. 

Deuteration is also strongly affected by changes in temper¬ 
ature. For example, the (late-time) deuteration degree of H.j is 
much lower at T - 20 K than at 7’ = 10 K because the higher 
temperature causes 1) the deuteration reactions (Ht + HD —» 
H 2 D + + H 2 etc.) to proceed more efficiently in the backward 
direction and 2) less severe depletion of other important reac¬ 
tion partners of H 2 D + such as CO and N 2 , observed also as en- 
chanced deuterations of HCO + and N 2 H + . Consequently, very 
little HD depletion is present at T - 20 K regardless of density. 

Although not readily apparent in the figures included here, 
we find a local minimum in deuteration at around 10 K if we 
study the abundances of deuterated species as a function of tem¬ 
perature but at fixed density. The temperature dependence is 
smoothed out if tunneling is included in the model, which also 
affects the deuteration degrees. This phenomenon might provide 
a means to investigate if tunneling is indeed effective on grain 
surfaces by studying the temperature dependence of the deu¬ 
terium fractionation. This will be the subject of a future paper. 

5. Conclusions 

We have developed new chemical reaction sets for gas-phase and 
grain-surface chemistry, using the osu_®l_20Q9 reaction set as 
a template. The new sets presented here are evolved versions 
of those presented in SI3; the present model includes deutera¬ 
tion of species with up to six atoms (so that we can include the 
formation pathways of, e.g., ammonia and methanol), and the 
spin-state chemistry for the species involved in the water and 
ammonia formation networks. 

We chose a simplified physical model, carrying out chemical 
calculations with a pseudo-time-dependent chemical code at dif¬ 
ferent densities and temperatures (but with otherwise fixed phys¬ 
ical parameters) in order to facilitate straightforward comparison 
of future modeling results against those presented in this paper. 
Calculations were performed with gas-phase and gas-grain mod¬ 
els. Special attention was given to the ortho/para ratios of water 
and ammonia, and to the effect of either including or excluding 
deuterium in the chemical model. 

We find that the o/p ratios of water and ammonia are ~ 1.6 
and ~ 0.5, respectively, for t > 10 5 yr irrespective of density or 
model (gas-phase or gas-grain) used. The o/p ratio of water is 
clearly lower than the value (~ 3) observed toward translucent 
clouds, but we find that this apparent disagreement disappears 
if we consider low values of visual extinction (as opposed to 
Am = 10 assumed elsewhere in the present paper). 

At early times, the o/p ratios vary with density, but little vari¬ 
ation is seen between gas-phase and gas-grain models, implying 
that the o/p ratios are determined by gas-phase processes. Our re¬ 
sults also show that, in the range of physical parameters consid¬ 
ered, excluding deuterium from the model has only a marginal 
effect on the abundances of non-deuterated species. 

We find that an increase of the temperature from 10 to 
20 K generally decreases the depletion of heavy molecules onto 
grain surfaces, and this decreases also the deuteration degree. 


However, nitrogen chemistry proceeds differently; ammonia and 
N 2 H + deplete more strongly at T = 20 K than at T = 10 K. 

The physical model of the present paper was deliberately 
kept simple in order to facilitate straightforward benchmarking 
of other models against the results presented here. However, we 
note that while the o/p ratio of ammonia did not show density de¬ 
pendence in our calculations, a more complete parameter-space 
study of ammonia (and water) spin-state chemistry, in the con¬ 
text of a gas-grain model, is called for. 
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Appendix A: On the late-time ortho/para ratios of 
ammonia and water 


In the interstellar medium, ammonia production is dominated 
by the electron recombination of the ammonium ion, NHJ, and 
its destruction occurs mainly through a charge transfer reaction 
with H + . At late times of chemical evolution, the proton transfer 
reaction with H j returns ammonia back to NHJ. The dominant 
reactions are shown schematically in Fig. lA. II The destruction 
rates of pNH3 and 0NH3 in reactions with H + , pH^, and 0H3 are 
equal, and consequently the 0/P-NH3 ratio is determined by the 
nuclear spin branching ratios in the electron recombinations of 
para-, ortho-, and meta-NHJ. One obtains the relationship 


0 /P-NH 3 = 


m/p-NH+ + l 0 /p-NH+ 
1 + fo/p-NHJ 


(A.l) 


Water is produced primarily in the electron recombination of 
the hydronium ion, H30 + , and like ammonia, destroyed mainly 
in the charge transfer reaction with H + . The recombination of 
pH30 + yields both 0H2O and pITO at equal ratios, whereas 
oH 3 0 + only yields ortho-water (along with hydroxyl). We ob¬ 
tain the relationship 

o/p-HiO = 2 o/p-H30 + + 1 . (A. 5 ) 


The production of H30 + is dominated by H20 + + H2 —» H30 + + 
H. The spin selection rules result in the relationship 


o/p-H 3 0 + 


|o/P-H 2 0 + 

1 + |o/p-H 2 0 + ' 


(A.6) 


(h 3 o + ) 



The ammonium ion is predominantly formed in the reactions 
pNH J + pfT and oNH j + pkb. These reactions determine the 
nuclear spin ratios of NH|, because the electron recombination 
rates are equal for the different nuclear spin species. The follow¬ 
ing relationships are obtained: 

m/p—NH 4 + = fo/p-NHJ 

o/p-NH+ = f o/p-NH++ f . ' ’ 

The primary production pathway to NH3, NH) + H2 —> 
N'H J + H, becomes at late times less important than the charge 
transfer reaction between NH3 and H + , and in this situation one 
obtains roughly equal o/p ratios for NH3 and NH3: 

o/p-NHj ~ 0/P-NH3 . (A. 3 ) 

Finally, substituting Eqs. (IA. 2 I) and (l/\. 3 b to Eq. ( IA. 1 I) . one ob¬ 
tains (with a little algebra) the steady-state ratio 

0/P-NH3 ~ 0.4 . (A. 4 ) 

This value is about 10 % lower than the one predicted by our 
simulation with the full reaction set. 

The ammonia abundance in interstellar molecular clouds is 
frequently derived using observations of the (1,1) and (2,2) in¬ 
version lines at A — 1.2 cm which both represent para-MT}. 
The total ammonia abundance is then derived by assuming 
0/P-NH3 = 1 or that the ortho and para states are populated 
according to LTE. The latter assumption implies at 10 K that 
o/p - NH3 = 3 . 3 . Our result suggests that these previous ob¬ 
servational estimates of the total ammonia abundance can be un¬ 
realistically large. 


Fig. A. 2 . Dominant reactions involving water and hydroxyl. 

The most important reactions determining the abundance 
of water and related molecules in the gas-phase are shown in 
Fig. lA . 21 The dominant pathway to FEO + is the H atom abstrac¬ 
tion reaction OH + + pFE —» (p or o)H20 + + H, which leads at 
equal probabilities to pH20 + and 0H2OA The second reaction 
of importance is the charge transfer reaction FEO + H + (with 
a 20% share of the H20 + production), and the third is proton 
transfer from H3 (~ 10 %). Omitting the secondary reactions, 
one obtains 

o/p-H 2 0 + ~ 1 . (A. 7 ) 

The full reaction set yields o/p-H20 + ~ 1.6 at late times. The 
substitution of Eqs. (IA.6I ) and (IA. 7 I) to Eq. (IA. 5 I ) yields a value 

o/p-H 2 0 ~ 1.4 , (A.8) 

representing steady state. This is about 20 % lower than the late¬ 
time o/p ratio in our simulation. This discrepancy can be under¬ 
stood by the fact that in the very simple model described above, 
we have omitted the production of FEO + from FEO and OH. 
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Appendix B: New branching ratios for water chemistry & updates to SI3 


Table B.l. The branching ratios of the most important reactions in the water formation network. The reactions given below are the 
most important in the sense that they are the dominant reactions for the various species, as confirmed by the output by our chemical 
code. The destruction of H 3 0 + by HCN and HNC is important at early times, but electron DR reactions dominate at late times. 




Chemical reaction 




Branching ratio 

oH+ 

0 

-> 

OH+ 

oH 2 


1.00 

PHI 

0 

—» 

OH+ 

pH 2 


0.50 

pii,; 

O 

-> 

OH+ 

oH 2 


0.50 

oH+ 

0 

-> 

oH 2 0 + 

H 


1.00 

pH; 

0 

— » 

pH 2 0 + 

H 


0.50 

pH; 

0 

-> 

oH 2 0 + 

H 


0.50 

OH+ 

oH 2 

— > 

pH 2 0 + 

H 


0.17 

OH+ 

oH 2 

- > 

oH 2 0 + 

H 


0.83 

OH+ 

pH 2 

—» 

pH 2 0 + 

H 


0.50 

OH+ 

pH 2 

- > 

oH 2 0 + 

H 


0.50 

oH 2 0 + 

oH 2 

-> 

oH 3 0 + 

H 


0.67 

oH 2 0 + 

oH 2 

— > 

pH 3 0 + 

H 


0.33 

oH 2 0 + 

pH 2 

-> 

oH 3 0 + 

H 


0.33 

oH 2 0 + 

pH 2 

—» 

pH 3 0 + 

H 


0.67 

pH,0 + 

oH 2 

-> 

oH 3 0 + 

H 


0.33 

pH 2 0 + 

oH 2 

— > 

pH 3 0 + 

H 


0.67 

pH 2 0 + 

pH 2 

—» 

pH 3 0 + 

H 


1.00 

oH+ 

oH 2 0 

-> 

pH 3 0 + 

pH 2 


0.07 

oH+ 

oH 2 0 

- > 

oH 3 0 + 

pH 2 


0.08 

oHJ 

oH 2 0 

—> 

oH 3 0 + 

oH 2 


0.62 

oHJ 

oH 2 0 

—> 

pH 3 0 + 

oH 2 


0.23 

0 H 3 

pH 2 0 

— » 

oH 3 0 + 

pH 2 


0.25 

oH+ 

pH 2 0 

- > 

oH 3 0 + 

oH 2 


0.25 

oH+ 

pH 2 0 

—» 

pH 3 0 + 

oH 2 


0.50 

pH; 

oH 2 0 

- > 

pH 3 0 + 

pH 2 


0.13 

pH; 

oH 2 0 

-> 

oH 3 0 + 

pH 2 


0.17 

pH; 

oH 2 0 

- > 

oH 3 0 + 

oH 2 


0.23 

pH; 

oH 2 0 

— > 

pH 3 0 + 

oH 2 


0.47 

pH; 

pH 2 0 

— > 

pH 3 0 + 

pH 2 


0.40 

pH; 

pH 2 0 

-> 

oH 3 0 + 

oH 2 


0.20 

pn; 

pH 2 0 

— » 

pH 3 0 + 

oH 2 


0.40 

oH 3 0 + 

e _ 

- > 

OH 

H 

H 

1.00 

pH 3 0 + 

e - 

—> 

OH 

H 

H 

1.00 

oH 3 0 + 

e _ 

- > 

oH 2 0 

H 


1.00 

pH 3 0 + 

e _ 

—» 

pH 2 0 

H 


0.50 

pH 3 0 + 

e _ 

-> 

oH 2 0 

H 


0.50 

oH 3 0 + 

e _ 

-> 

OH 

oH 2 


1.00 

pH 3 0 + 

e - 

—» 

OH 

pH 2 


0.50 

pH 3 0 + 

e _ 

-> 

OH 

oH 2 


0.50 

oH 3 0 + 

e _ 

-> 

oH 2 

H 

0 

1.00 

pH 3 0 + 

e _ 

—> 

pH 2 

H 

0 

0.50 

pH,0 + 

e _ 

-> 

oH 2 

H 

O 

0.50 

oH 3 0 + 

HNC 

-> 

HCNH + 

pH 2 0 


0.12 

oH 3 0 + 

HNC 

-> 

HCNH+ 

oH 2 0 


0.88 

pH 3 0 + 

HNC 

— > 

HCNH + 

pH 2 0 


0.38 

pH 3 0 + 

HNC 

-> 

HCNH+ 

oH 2 0 


0.62 

oH 3 0 + 

HCN 

— > 

HCNH + 

pH 2 0 


0.12 

oH 3 0 + 

HCN 

-> 

HCNH+ 

oH 2 0 


0.88 

pH 3 0 + 

HCN 

—> 

HCNH + 

pH 2 0 


0.38 

pH 3 0 + 

HCN 

—» 

HCNH+ 

oHiO 


0.62 
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Table B.2. Rate coefficient data from lLe Gal et al.l (2014) included in the present model. We have included some of the Le Gal et 
al. data “as is”, while for most reactions we have adopted a total rate coefficient in which contributions from different spin states 
are added together. When these reactions are separated with our spin-state separation routine, we recover the coefficients given in 
iLeGal et al.l d20l4 . 




Chemical reaction 




a 

P 

r 

As is 

N + 

oH 2 

-> 

NH + 

H 


4.20 x 10 10 

-0.15 

44.10 

N + 

pH 2 

— > 

NH + 

H 


8.35 x 10 - 10 

0.00 

168.50 

H + 

oHi 

— > 

H + 

pH 2 


1.82 x 10~ 10 

0.13 

- 0.02 

H + 

pH 2 

— » 

H + 

oH 2 


1.64 x 10 ~ 9 

0.13 

170.50 

HCO+ 

oH, 

— > 

HCO+ 

pH 2 


1.27 x 10 10 

0.00 

0.00 

HCO + 

pH 2 

— » 

HCO + 

oH 2 


1.14 x 10 ~ 9 

0.00 

170.50 

N 

OH 

— » 

NO 

H 


8.90 x 10-“ 

0.20 

0.00 

N 

NO 

— > 

n 2 

0 


7.20 x 10-“ 

0.44 

12.70 

N 

CN 

— » 

n 2 

C 


8.80 x 10” 11 

0.42 

0.00 

N 

CH 

— > 

CN 

H 


1.70 x 10 10 

0.18 

0.00 

C 

NO 

— > 

CN 

0 


6.00 x 10 -“ 

-0.16 

0.00 

C 

NO 

—> 

CO 

N 


9.00 x 10 ” 11 

-0.16 

0.00 

Total 

NH + 

h 2 

- > 

NHf 

H 


1.28 x 10~ 9 

0.00 

0.00 

NH( 

h 2 

—» 

NHt 

H 


2.70 x 10 ~ 10 

0.00 

0.00 

NHf 

h 2 

—> 

NH+ 

H 


2.40 x 10 - 12 

0.00 

0.00 

HJ 

NH 

—> 

NH( 

h 2 


1.30 x 10 ~ 9 

0.00 

0.00 

H : t 

nh 2 

—» 

NH+ 

h 2 


1.80 x 10~ 9 

0.00 

0.00 

H;t 

nh 3 

—> 

NH+ 

h 2 


9.12 x 10 ~ 9 

0.00 

0.00 

HCO+ 

NH 

—> 

NH( 

CO 


6.40 x 10 10 

0.00 

0.00 

HCO+ 

nh 2 

—> 

NH+ 

CO 


8.90 x 10 10 

0.00 

0.00 

HCO + 

nh 3 

—> 

NH+ 

CO 


1.92 x 10 9 

0.00 

0.00 

N 2 H + 

NH, 

—> 

NH+ 

n 2 


2.30 x 10 ~ 9 

0.00 

0.00 

n 2 h + 

e - 

—> 

n 2 

H 


2.77 x 10 ~ 7 

-0.50 

0.00 

n 2 h + 

e _ 

—> 

NH 

N 


2.09 x 10 ~ 8 

-0.50 

0.00 

NH( 

e _ 

—> 

NH 

H 


1.17 x 10 ~ 7 

-0.50 

0.00 

NH( 

e _ 

—» 

N 

H 

H 

1.71 x 10~ 7 

-0.50 

0.00 

NH( 

e _ 

—> 

N 

h 2 


1.20 x 10~ 8 

-0.50 

0.00 

NH} 

e~ 

—> 

nh 2 

H 


1.55 x 10 ~ 7 

-0.50 

0.00 

NH+ 

e _ 

—» 

NH 

H 

H 

1.55 x 10 ~ 7 

-0.50 

0.00 

NHJ 

e _ 

—> 

NH, 

H 

H 

1.22 x 10- 7 

-0.60 

0.00 

NH+ 

e~ 

—» 

nh 2 

h 2 


1.88 x 10“ 8 

-0.60 

0.00 

NH) 

e _ 

—» 

NH, 

H 


8.00 x 10- 7 

-0.60 

0.00 

H 3 0 + 

e - 

—> 

OH 

h 2 


3.00 x 10 ~ 8 

-0.50 

0.00 

h 3 o + 

e~ 

—» 

OH 

H 

H 

2.60 x 10~ 7 

-0.50 

0.00 

h 3 o + 

e _ 

—> 

H,0 

H 


1.10 x 10~ 7 

-0.50 

0.00 

h 3 o + 

e _ 

—> 

h 2 

H 

0 

2.80 x 10~ 9 

-0.50 

0.00 
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Appendix C: Calculations at different temperatures 

In this appendix, we present the results of calculations otherwise similar to those presented in Sect. [3] but produced assuming either 
^gas = 7du St = 15 K or r gas = r dus t = 20K. 
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Fig. C.l. As Fig. [2] but calculated with 7 gils = 7’ ( | ust = 15 K. 
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Fig. C.2. As Fig. |3] but calculated with 7' g;ls = rd ust = 15 K. 
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Fig. C.5. As Fig. |3] but calculated with 7 g;ls = 7’ ( | ust = 20 K. 
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Fig. C.6. As Fig. [4] but calculated with 7 g;ls = rd us , = 20 K. 
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Appendix D: Calculations without deuterium 

The inclusion of deuterium in a chemical model may decrease the abundances of nondeuterated species especially at high densities 
in the presence of depletion, where deuterium fractionation is strong. To investigate how our results change when deuterium is 
excluded, we have run model calculations at T = 10 K setting the initial HD abundance (i.e., the initial deuterium abundance) to 
zero. The results of these calculations are presented in Figs. iDTltolDJl which should be compared against Figs.[2]to[4]in the main 
text. 

It is observed that in the gas-phase model, where there is little deuteration regardless of density, exluding deuterium from the 
calculations makes virtually no difference in the abundances of the non-deuterated species. Even in the gas-grain model, only very 
small enhancements are observed in the abundances of non-deuterated species at high density. Also, the o/p ratios of HiO and NH 3 
are only very slightly modified in the gas-grain model, and practically not at all in the gas-phase model. Our models thus imply that 
disregarding deuteration in chemical models when T > 10 K will likely not lead to large errors in the abundances of non-deuterated 
species. We note that we have not carried out a full parameter-space exploration of this issue here, and have only explored the effect 
of density (and temperature; Appendix^ on the results. 



Fig.D.l. Abundances of selected species in a gas-phase model excluding deuterium. This figure should be compared against Fig. [2] 



Fig.D.2. Abundances of selected species in a gas-grain model excluding deuterium. This figure should be compared against Fig. [3] 
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Fig.D.3. As Fig. [4] but excluding deuterium. 
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Appendix E: Calculations with quantum tunneling included 

Below, we present the results of calculations performed at T — 10 K, but including quantum tunneling on grain surfaces. When 
tunneling is included, the thermal diffusion rate defined by Eq. ((hj) is replaced by the tunneling diffusion rate dHasegawa et al.lfl992h 

< ff ' q = exp [-2 (a/fi) (2 m k B E? s ) l/2 ] , (E. 1) 

where a is the width of the (rectangular) tunneling barrier. We assume a = 1 A. Also, the reaction probability k,j is in the presence 
of tunneling replaced by 

^ = exp [-2 (a/ti) (2fik B E a ) 1 ' 2 ] , (E.2) 

where /./ is the reduced mass of the reactants (not to be confused with the mean molecular weight of the gas as defined in Eq.0. The 
reaction rate coefficient assumes the same form as when tunneling is excluded (Eq.[8j. In the present model, we only allow atomic 
H and D to tunnel, i.e., Eqs. (IE. II) and (IE.21) are only used for those reactions where either of these species is present as a reactant. 

Figure IETI presents the results of calculations at 7’ = 10 K with quantum tunneling included. This figure should be compared 
against Fig. [2 in the main text. Evidently, tunneling influences our results only at long timescales (> 10 6 yr), and even then the 
influence on deuteration for example is small. The differences between the models arise because tunneling allows reactions with 
activation barriers to proceed efficiently. We refer the reader to S13 for more discussion on tunneling and its effects on deuteration. 
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Fig.E.l. As Fig. [3] but with quantum tunneling included. 
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